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ABSTRACT 

Context. The physics of the pulsar magnetosphere near the neutron star surface remains poorly constrained by observations. Although 
about 2000 pulsars have been discovered to date, little is known about their emission mechanism, from radio to high-energy X-ray 
and gamma-rays. Large vacuum gaps probably exist in the magnetosphere, and a non-neutral plasma partially fills the neutron star 
surroundings to form an electrosphere. 

Aims. In several previous works, we showed that the differentially rotating equatorial disk in the pulsar's electrosphere is diocotron 
unstable and that it tends to stabilise when relativistic effects are included. However, when approaching the light cylinder, particle 
inertia becomes significant and the electric drift approximation is violated. In this paper, we study the most general instability, i.e. by 
including particle inertia effects, as well as relativistic motions. Electromagnetic perturbations are described in a fully self-consistent 
manner by solving the cold-fluid and Maxwell equations. This general non-neutral plasma instability is called the magnetron instability 
by plasma physicists. 

Methods. We linearise the coupled relativistic cold-fluid and Maxwell equations. The non-linear eigenvalue problem for the perturbed 
azimuthal electric field component is solved numerically with standard techniques for boundary-value problems like the shooting 
method. The spectrum of the magnetron instability in a non-neutral plasma column confined between two cylindrically conducting 
walls is computed for several cylindrical configurations. For a pulsar electrosphere, no outer wall exists. In this case, we allow for 
electromagnetic wave emission propagating to infinity. 

Results. First we checked our algorithm in the low-density limit. We recover the results of the relativistic diocotron instability. When 
the self-field induced by the plasma becomes significant, it can first increase the growth rate of the magnetron instability. However, 
equilibrium solutions are only possible when the self-electric field, measured by the parameter and tending to disrupt the plasma 
configuration, is bounded to an upper limit, ^e niax- For .Se close to but smaller than this value s^ y„.^^, the instability becomes weaker or 
can be suppressed as was the case in the diocotron regime. 

Conclusions. When approaching the light-cylinder, particle inertia becomes significant in the equatorial disk of the electrosphere. 
Indeed, the rest-mass energy density of the plasma becomes comparable to the magnetic energy density. The magnetron instability 
sets in and takes over the destabilisation of the stationary flow initiated by the diocotron instability close to the neutron star surface. 
As a consequence, the flow in the pulsar inner magnetosphere is highly unstable, leading to particle diffusion across the magnetic field 
line. Therefore, an electric current can circulate in the closed magnetosphere and feed the wind with charged particles. 

Key words. Instabilities - Plasmas - Methods: analytical - Methods: numerical - pulsars: general 



1. INTRODUCTION 

This year, we celebrate the 40th year of the discovery of the first 
pulsar. Nevertheless, the detailed structure of charge distribution 
and electric-current circulation in the closed magnetosphere of a 
pulsar remains poorly understood. Although it is often assumed 
that the plasma fills the space entirely and corotates with the neu- 
tron star, it is on the contrary very likely that it only partly fills 
it, leaving large vacuum gaps between plasma-filled regions. The 
existence of such gaps in aligned rotators has been very clearly 
established by Krause-Polstorff & Michel ( 1985a b). Since then, 
a number of different numerical approaches to the proble m have 
confirmed their conclusions, including some work by Rvlovl 
(1989), Shibata (1989), Zachariades (1993), Neukirchf (Il993l) 



fThielheim & Wolfsteller (1994), iSpitkovskv & AronsI (I2OO2I) . 
and ourselves (Petri et al., 2002bi This conclusion about the ex- 
istence of vacuum gaps has been reached from a self-consistent 
solution of the M axwell equat i ons in the case of the aligned 
rotator. Moreover, ISmith et al] (1200 lb have shown by numer- 
ical modelling that an initially filled magnetosphere like the 
Goldreich-Julian model evolves by opening up large gaps and 
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Stabilise s to the partially filled and part ially void sol ution found 
by Kra use-Polstorff & Michell ( Il985ah and also by iPetri et alj 
(i2002bj) . The status of models of the pulsar magnetospheres, or 
electro spheres, has recently been critically reviewed by MicheJ 
(l2005h . A solution with vacuum gaps has the peculiar property 
that those parts of the magnetosphere that are separated from the 
star's surface by a vacuum region are not corotating and so suf- 
fer differential rotation, an essential ingredient that will lead to 
non-neutral plasma instabilities in the closed magnetosphere, a 
process never addressed in detail. 

This raises the question of the stability of such a charged 
plasma flow in the pulsar magnetosphere. The differential ro- 
tation in the equatorial, non-neutral disk induces a non- neutral 
plasma instability that is well known to plasma physicists dOneill 
1980; Davidson 1990; O'Neil & Smith 1992). Their good con- 
finement properties (trapped particles can remain on an almost 
unperturbed trajectory for thousands of gyro-periods) makes 
them a valuable tool for studying plasmas in laboratory, by using 
for instance Penning traps. In the magnetosphere of a pulsar, far 
from the light cylinder and close to the neutron star surface, the 
instability reduces to its non-relativistic and electrostatic form, 
the diocotron instability. The linear development of this insta- 
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bility for a differe ntially rotating charged disk was studied by 
iPetri et al.1 ( l2002ab . in the thin disk Umit, and by 'Petril (l2007albt) 
in the thick disk limit. It both cases, the instability proceeds at a 
growth rate comparable to the star's rotation ra te. The non-linea r 
development of this instability was studied bv .Fetri et al.l(l2003l) . 
in the framework of an infinitely thin disk model. They have 
shown that the instability causes a cross-field transport of these 
charges in the e quatorial disk, evolving into a net out-flowing 
flux of charges. ISpitkovskv & AronsI (|2002|) have numerically 
studied the problem and concluded that this charge transport 
tends to fill the gaps with plasma. The appearance of a cross- 
field electric current a s a result o f the diocotron instability has 
been observed by Pasq uini & Faians (200 2) in laboratory exper- 
iments in which charged particles were continuously injected in 
the plasma column trapped in a Malmberg-Penning configura- 
tion. 

A general overview of the equilibrium and stability proper- 
ties of non-neutral p lasmas in Carte s ian an d cylindrical geome- 
try can be found in [Davidson et al.l (Il991h . iTsang & DavidsonI 
( 1198 6') describe how to compute fully self-consistent general 
equilibria configuration for a cold-fluid plasma in a cylindri- 
cal diode. This is useful to investig ate the stability propertie s 
in magnetron devices as presented in Davidson & Tsangl(ll986l) . 

The aim of this work is to extend the previous work bv lPetril 
( l200 7a"b) on the diocotron instability by including particle in- 
ertia effects. In this paper we present a numerical analysis of 
the linear growth of the relativistic magnetron instability for a 
non-neutral plasma column. The paper is organised as follows. 
In Sect. 12] we describe the initial setup of the plasma column 
consisting of an axially symmetric equilibrium between two con- 
ducting walls. We give several equilibrium profiles useful for the 
study of the magnetron instability in different configurations. In 
Sect. [3] the non-linear eigenvalue problem satisfied by the per- 
turbed azimuthal electric field component is derived. The algo- 
rithm to solve the eigenvalue problem is checked against known 
analytical results in the low-density limit (diocotron instability). 
Sect, m Then, applications to some typical equilibrium configu- 
rations are shown in Sect. |5] First we consider a plasma column 
with constant density. Next, we study the effect of the cylindrical 
geometry (curvature of the flow) and the transition to the planar 
diode limit. Finally, the consequences of the magnetron instabil- 
ity on the pulsar electrosphere is investigated. The conclusions 
and the possible generalisation are presented in Sect. |6] 



2. THE MODEL 

We study the motion of a non-neutral plasma column of infinite 
axial extend along the z-axis. We adopt cylindrical coordinates 
denoted by (r, {p, z) and define the corresponding orthonormal 
basis vectors by (ei-,e^,ez)- The geometric configuration is the 
same as the one described in our previous works (Petri 2007 a. 3). 

In this section, we briefly summarise the equilibrium con- 
ditions imposed on the plasma and give some typical examples 
of equilibrium configurations for specified velocity, density and 
electric field profiles. 

We consider a single-species non-neutral plasma consist- 
ing of particles with mass me and charge q trapped between 
two cylindrically conducting walls located at the radii Wi and 
W2 > Wi. The plasma column itself is confined between the radii 
R\ > Wi and R2 < W2, with Ri < R2- This allows us to take into 
account vacuum regions between the plasma and the conducting 
walls. Because we solve the full set of Maxwell equations, there 
is also the possibility for the plasma to radiate electromagnetic 



waves to infinity. In order to take this effect into account, we re- 
move the outer wall if necessary and replace it by outgoing wave 
boundary conditions at R2. 

2. 1 . Equations of motion 

Let us describe the equation of motion satisfied by the plasma 
column. Each particle evolves in the self-consistent electromag- 
netic field partly imposed by an external device and partly in- 
duced by the plasma itself. The motion of the plasma column 
is governed by the cold-fluid equations, i.e. the conservation of 
electric charge and the conservation of momentum, respectively. 



dt 



-Hdiv(peV) = 



(1) 



-+v-—\{ym,v) = q(E + vAB) (2) 

We introduced the usual notation, pe for the electric charge den- 
sity, V for the velocity, (E, B) for the electromagnetic field. The 

Lorentz factor of the flow is 7 = 1 / -y/l'^-'vVc^- Eqs. ([T]i and (|2]l 
are supplemented by the full set of Maxwell equations, 

roiE = - — 

dt 



rotB = Ha] + bqiiq — 
ot 

divE = ^ 

£() 

divB = 



(3) 

(4) 

(5) 
(6) 



where yUo is the magnetic permeability and sq the electric permit- 
tivity. For a non-neutral plasma, the current density corresponds 
to convective motion in the plasma and is therefore related to the 
charge density by 

J = Pe V (7) 

The set of Eqs. ([B-© represents the most general treatment 
of the cold fluid behavior, including relativistic and electromag- 
netic effects as well as inertia of the particles (the electric drift 
approximation, used to study the diocotron instability, is re- 
placed by Eq. (|2|). 

2.2. Equiiibrium of the plasma column 

In equilibrium, the particle number density is ne('") and the 
charge density is Ps{r) = q n^{r). Particles evolve in a cross elec- 
tric and magnetic field such that the equilibrium magnetic field 
is directed along the z-axis whereas the equilibrium electric field 
is directed along the r-axis. On average, the stationary motion is 
only azimuthal. The electric field induced by the plasma itself is 



(8) 



The magnetic field is made of two parts, an imposed external 
applied field, B^, assumed to be uniform in the region outside 
the plasma column, and a plasma induced field. Bp such that the 
total magnetic field is 

B = Bp-HBo = Bzez. (9) 

Therefore, for azimuthally symmetric equilibria, the steady-state 
Maxwell-Gauss and Maxwell- Ampere equations satisfy 

r or £() 



SB, 
dr 



-y"OPe VV- 



(10) 
(11) 
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In the stationary state, d/dt -0, the balance between the Lorentz 
force and the centrifugal force for a fluid element in Eq. (|2| is 
expressed as 



7 me — +q (E, + B^) - 0. 



(12) 



It is convenient to introduce the non-relativistic plasma and cy- 
clotron frequencies respectively by 



Peg 

nig So 
nig 



The corresponding relativistic expressions are respectively 



2 P 

W = 

P y 



(13) 
(14) 

(15) 
(16) 



axial magnetic field, Bz, the charge density, pe, and the azimuthal 
speed of the guiding center, v^. Prescribing one of these profiles, 
the remaining three are found self-consistently by solving the set 
of Eqs. ( [Tol l. (fTTI) and (fTSI) . In the next subsections, we show how 
to derive these quantities for some typical examples in which ei- 
ther the velocity profile, the density profile or the electric field 
are imposed. 



2.3. Specified velocity profile 

Let us first assume that the velocity profile = rQ is pre- 
scribed, Q being the angular velocity at equilibrium. This case is 
well-suited for the study of the pulsar's electrosphere in which 
the plasma is in differential rotation and evolves in a dipolar 
magnetic field. As already noticed before, the differential ro- 
tation is essential to the presence of the instability. The other 
equilibrium quantities, (Er,B^,ps), are easily derived from Q. 
Indeed, inserting pe from Eq. (fTOl l and Ei from Eq. ( fT2l l into 
Maxwell-Ampere equation ( fTTI ). the magnetic field satisfies a 
first order linear ordinary diff'erential equation 



where y = 1 / ^1 - v^/c^ corresponds to the Lorentz factor of 

the azimuthal flow. The plasma equilibrium is governed by two 
antagonistic effects. On one hand, the electric field induced by 
the plasma itself exerts a repelling force on the fluid element, 
trying to inflate the plasma. On the other hand, the magnetic 
field confines the plasma by incurving the particle trajectories. 
To quantify the strength of each effect, focusing magnetic field 
and defocusing electric field, it is useful to introduce the self- 
field parameter defined by 



(17) 



Note that this self-field parameter is a local quantity defined on 
every point in space as are the plasma and cyclotron frequencies. 
It can therefore depend on the radius r and is generally not a 
constant throughout the plasma column. We should write it Se(r) 
but in order to avoid overloading, we do not explicitly show the 
space dependence of any quantity as they should all depend on 
r except otherwise specified. The diocotron regime, or electric 
drift approximation, corresponds to the low-density limit, i.e. to 
ie ^ 1 in the whole plasma column. We assume that the electric 
field induced by the plasma vanishes at the inner wall, at r = Wi, 
i.e. 



Integrating Eq. (fTOl i therefore gives for the electric field gener- 
ated by the plasma. 



EJr)^— I p,{r')r' dr' e. 



Eor Jw, 



(19) 



For the magnetic field induced by the plasma, we solve Eq. (fTTT i 
with the boundary condition B^iRi) - Bq. This simply states that 
the total magnetic field outside the plasma column has to match 
the magnetic field imposed by an external device. Integrating 
within the plasma column, it is given by : 



B^(r) = Bo - jiQ 



Jr. 



(20) 



Any equilibrium state is completely determined by the fol- 
lowing four quantities, the total radial electric field, E^, the total 



dB, 
dr 



r 



d me c d / ^\ 

or q or ^ ' 



The Lorentz factor of the flow is 
1 



r = 



rQ 

= — 

c 



(21) 

(22) 
(23) 



From Maxwell-Ampere equation, Eq. ( fTTT l. the charge density is 
found by 



Pe = - 



1 dB, 



fiQV^ dr 



(24) 



The electric field is recovered from the force balance equation, 
Eq. O, 



ywe v; 

Er = Bz 

qr 



(25) 



(18) 2.4. Specified density profile 



Another interesting case corresponds to a specified charge den- 
sity profile Pe. Replacing Eq. ( fT9l ) and d20] i into the force balance 
Eq. ( fTSI i. the rotation profile is solution of a non-linear Volterra 
integral equation 



Vl -r2Q2/c2 



if 

f Jw, 



Q^{r')r'dr' 



(26) 



+ rO 



-4 r al{r')D.{r')r'dr' 
c Jr2 



= 



Qco = qBo/nifi is the non-relativistic cyclotron frequency asso- 
ciated to the external magnetic field Bq. Knowing Q, the same 
procedure as in the previous subsection for a specified velocity 
profile is applied, i.e. the electromagnetic field is calculated ac- 
cording to Eq. ( l2TT i and ( |25] l. 
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2.5. Specified electric field 

It is also possible to specify the equilibrium radial electric field. 
An interesting case is given by 



Er(r) = 







-c B, 



sinhci'(r - Ri) 

cosha(R2-Rl) 
Ri IwhaiRo- R\) 



Wx<r<Rx 
Ri <r<R2 

R2 < r <W2 



(27) 



a is a constant useful to adjust the maximal speed of the column. 
The equilibrium electric field profile, Eq. dZTl i. enables us to in- 
vestigate the influence of the cylindrical geometry compared to 
the planar diode geometry. Indeed, in the limit of small curva- 
ture, R2 - R\ « ^1, the eigenvalue problem in cylindrical ge- 
ometry reduces to the planar diode case. The magnetic field is 
solution of an ordinary differential equation 



or rc or 



(28) 



The density is easily found from Maxwell-Poisson Eq. ( fTOb to 
be 



noq 



(29) 



£() Bo a 



cosha(R2-Rl) 



cosh Of (r - RI) + 



sinhQ'(r - RI) 



where no is the plasma density within the column. Finally, the 
velocity is given by the centrifugal and Lorentz forces balance, 
Eq. ([nil- 

3. LINEAR ANALYSIS 

In this section, we derive the eigenvalue problem for the mag- 
netron instability in the most general case. We apply the standard 
linear perturbation theory. All scalar perturbations of a physical 
quantity X like electric potential, density, and velocity compo- 
nents, are expressed by the expansion 

X(r,if,t)^X(r)e'^'^-"'^ (30) 

where I is the azimuthal mode and co the eigenfrequency. A more 
general treatment of the perturbation in full 3D allowing for a 
stratified vertical structure of the plasma column would need to 
introduce a wavenumber k, directed along the z-axis, in the phase 
term, like (kz + lip - cuit). Techniques to deal with this general 
expansion are discussed in Sect. 3.1 o f lPetril(l2007al ). 

The decomposition in Eq. ( |30] l only allows for a global mode 
to propagate in the plasma. The frequency does not depend on 
the radius. Nevertheless, we could imagine to split the column of 
plasma in diff'erent annular layers L, , let us say layers such that 
i G [[1--A^]]. Each of them possesses its own eigenfrequency w,. 
Then, by using matching conditions at the interface between suc- 
cessive layers, we could solve the eigenvalue problem and com- 
pute the radius dependent growth rates. This would be a general- 
isation of the technique employed to match the vacuum solution 
to the plasma column as described below. 



3.1. Linearisation of Maxwell equations 

We study the stability of the plasma column around the equi- 
librium mentioned in Sect. |2] An expansion to first order of the 
electromagnetic field around this equilibrium (£", B^) leads to 



E = E° + 6E 

= B? + 5B. 



(31) 
(32) 



and the same for the charge and current densities 
j = / + Sj (33) 

Pe = P° + 5pe (34) 

Linearising the set of Maxwell equations, ©-([Sll, we have 
eo 

/ <jj 6Bj 



1 d I 

- —{r6E,) + I - 6E^ 

r or r 



1 d I 
- —(r5E^)-i-5E, 
r or r 

i - 6B, 
r 

d 

-^5B, 
or 



= yuo Sj, - i 5E, 



The current density perturbation is 

6jr = Pe<5v,- 



(35) 

(36) 
(37) 
(38) 



(39) 
(40) 



It is convenient to introduce a potential 5E^p - -i I (pj r (we em- 
phasise that this function is not the scalar potential from which 
the electric field could be derived, but is just a convenient auxil- 
iary variable) such that the electric and magnetic field become 



6E, 



d4> r^ 
or P 



(jjr 

6B^ = — K 
I 



1 d(f) .juo 
c or a> 



1 



l-(ojr/lc)^ 
Maxwell-Gauss equation, is therefore written 

I d I d4>\ f 6p^ u> d Ir^ 



r dr 



dr 



(41) 
(42) 
(43) 

(44) 



In order to find the eigenvalue equation satisfied by the func- 
tion 0, we have to relate the density and velocity perturbations, 
6p^, 5v, and 5v^, to (j). These expressions are derived in the next 
paragraph. 

3.2. Linearisation of the fluid equations 

Linearising the conservation of momentum, Eq. (|2]), the radial 
and azimuthal perturbations in velocity are related to the pertur- 
bation in electromagnetic field as follows 

— i(Lo — ID.) 6\\ — (a>c + 2 Qb) 6v,p = 

{SE, + rQ.SB^)-i{cD-lSl)y^6v^ 



1 d 



-(r'ymSv, 



■6E^ 



(45) 



y r dr / ymg. 

Following lDavidson & Tsarigl(ll986h . we introduced the Coriolis 
frequency 



1 +y2 

From Eqs. dTIT l and (l42T l. we find 

5E,+rQ.6B^ = - 



(46) 



U)r Q.r\ d(f> r 

1 - — -iyU0Pe-o-(w-m)^Vr 

Ic c or 



(47) 
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From the continuity equation, ([T]), we get 
1 



(5pe = 



1 d I 

r or r 



iiuj-lQ.) 
Introducing the function 

/ \ d 
A = (We+2Qb) Wc + — ^(rV^^) 



(48) 



y r dr 



1+K 



2 2 



(49) 



we solve for the velocity perturbation as 
iq 



5v 



y Me A 



I o) r Q. r\ dd) 



- (wc + 2nb)-0 

r 



(50) 



1 



y me A 
- (w-/Q) 



^ 1 



cor Q.r 
Ic c 

, .2 a\ 



]_ d_ 
y r dr 



(r^Q) 



d4> 
Yr 



1 +K 



P c2 



(51) 



Therefore, the velocity perturbations are expressed in terms of 
the function (p as well as the density perturbation, Eq. ( l48T l. 

The eigenvalue problem for (f> is obtained by inserting 
Eqs. (l48b . dSOl l and ( l5TT l into Eq. (l44l) . After some algebraic ma- 
nipulations, we arrive at the final expression 



1 d_ 
r dr 



mil 

or 



Id) I Lor Q.r 
K 11 

(w-/0)r \ Ic c 



d_ 

dr 



P 

- — (1 +X.fi)4> 



-^(t^e+2Qb) 

A 



with the two auxiliary functions 



Xr 



2 2 



1 - 



a> r Q.r 
Ic c 



1 



Lo^r 



P, 



+ 2 (wc + 2 Qb) 



(52) 



(53) 



(54) 



The eigenvalue equation (l52T i is very general dPavidson & Tsan^ 
11986) . It describes the motions of small electromagnetic pertur- 
bations around the given equilibrium state, Eqs. (fTOl i and ( fTTl i. 
when inertia is taken into account. Many aspect of the relativis- 
tic magnetron instability can be investigated with this eigenvalue 
equation. In order to solve the eigenvalue problem, boundary 
conditions need to be imposed at the plasma/vacuum interface. 
They play a decisive role in the presence or absence of the insta- 
bility. How to treat the transition between vacuum and plasma or 
outgoing wave solutions is discussed in the next subsection. 

3.3. Boundary conditions 

In laboratory experiments, the plasma is usually confined be- 
tween inner and outer conducting walls. However, in pulsar elec- 
trospheres, no such outer device exists to constrain the electric 
field at the outer boundary. Radiation from the plasma could 
propagate into vacuum to infinity, carrying energy away from 
the plasma by Poynting flux. To allow for this electromagnetic 
wave production by the instabilities studied in this work, the 



outer wall is removed. The electromagnetic field is solved an- 
alytically in vacuum and matched to the solution in the plasma 
at the plasma/vacuum interface located at r - R2. First we dis- 
cuss the situation in which an outer wall exists and then consider 
outgoing waves. 

3.3.1. Outer wall 

When vacuum regions exist between the plasma column and the 
walls, special care is required at the sharp plasma/vacuum transi- 
tions. Indeed, the right-hand side of Eq. i5% then involves Dirac 
distributions 6{r) because the function f(r) = oj^ + 2 Ob)/ A 
is discontinuous at the edges of the plasma column, at Ri and 
i?2- Moreover, outside the plasma column, this function vanishes 
such that f(r) - for r < Ri and r > R2. 

In other words, its derivative has to be computed as 



— = 6{r-Ri) — 
or or 



-5{r-R2) % 

r=R, or 



df 
dr 



regular 



(55) 



where Ireguiar means the regular (or continuous) part of the deriva- 
tive, i.e. which does not involve 6 distributions. It vanishes in the 
vacuum regions, r < Ri and r > R2. Therefore, the first order 
derivative of (p is not continuous at these interfaces. To over- 
come this difficulty, we decompose the space between the two 
walls into three distinct regions: 

- region I: vacuum space between inner wall and inner bound- 
ary of the plasma column, with the solution for the electric 
potential denoted by (pi, defined for Wi < r < Ri ; 

- region II: the plasma column itself located between Ri and 
R2, solution denoted by (pu, defined for Ri < r < R2 ; 

- region III: vacuum space between the outer boundary of the 
plasma column and the outer wall, solution denoted by (pm, 
defined for R2 < r < W2. 

In regions I and III, the vacuum solutions satisfy the required 
boundary conditions, <pi{Wi) - and <pm{W2) = 0. 

The jump in the first order derivative d(p/dr at each 
plasma/vacuum interface are easily found by multiplying 
Eq. ( |52] i by r and then integrating around each discontinuity. 
Performing the calculation, we find the jump at /?i to be 



dr dr a>-lL}.{Ri) 



^ Q(Ri)RiCjRA 
1 : — X 



Ic 



A(Ri) Riil+XriRi)) 
Similarly, at the outer interface at R2, we obtain. 



(56) 



dr ^ ' dr^ ' cj-miR2) 



Ic 



^(^2) co,(R2) + 2Q.b(R2) 

A{R2) R2(l+Xr(R2)) 



(57) 



3.3.2. Outgoing wave solution 

Because of the wall located at r = W2, the outer boundary con- 
dition <pui{W2) - enforces E^{W2) - 0. It therefore prevents 
waves escaping from the system due to the vanishing outgoing 
Poynting flux, B./jiq = 0. In pulsar magnetospheres, no such 
wall exists. Thus, in order to let the system produce outgoing 
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electromagnetic waves, we remove the outer wall in this case 
and solve the vacuum wave equation for 0, which then reads 

^ ■ -cf>^0. (58) 



1 

r dr 



rK{r,(jS) — 
or 



This equation can also be derived directly from the vector wave 
equation 



(59) 



projected along the axis. To find the right outgoing wave 
boundary conditions, it is therefore necessary to solve the vector- 
wave equation in cylindrical coordinates usingjyector cylindri- 
cal harmonics as d e scribe d for instance in IStrattod (fl94lh and 
orse & FeshbachI (Il953l) . The solutions for the function (p to 
be an outgoing wave in a vacuum outside the plasma column, 
which vanishes at infinity, is given by (region III with W2 = +00) 

where the cylin drical outgoing wave function is given by Hi 
(IStrattonL[T94l . the Hankel function of the first kind and of or- 
der / related to the Bes s el f unctions by Hi(x) = Ji(x) + i Yi(x), 
jAbramowitz & Stegunl[T96 5). The prime ' means derivative of 
the function evaluated at the point given in parentheses, and K 
is a constant to be determined from the boundary condition at 
R2. Eliminating the constant K, we conclude that the boundary 
condition to impose on (p is 



L0R2 ^ 



R 



m4] 



liRi) - 



fill 



dr 



(Ri) = 0. 



(61) 



The boundary conditions expressed in region II for (pu are found 
by replacing 0jy(7?2) from Eq. (|57] |. Since <p is continuous, 
MM = (piM = (f>(R2)- We find 
f(jjR2\ UR2 „/a>R2 
c 



R2 



c 
d(pii 
dr 



<P(R2) ■ 



(R2) 



l(f>(R2) 

oj-m(R2) 



X 



X 1- 



Q.(R2)R2 CJR2\ '^l(R2) Wc(/?2) + 2 Qb(/?2) 



Ic I AiR2) R2{1 +XAR2)) 



= 0. (62) 



3.4. Set of annular layers 

Because the instability is related to strong velocity gradients, it 
is expected that the instability will only exist in the vicinity of 
this shear. In this paragraph, we briefly describe how to compute 
the radius dependent growth rate without going into details. 

As was already done for the plasma/vacuum interface, the 
plasma column is divided into a set of annular layers L,, each 
of them having their own eigenfrequency w, and a radial exten- 
sion from 7? J to > ^1 ■ Because no surface charge accumu- 
lates on the interfaces between successive layers, the matching 
conditions require the unknown function (p to be continuous as 
well as its first derivative when crossing the interface. Moreover, 
to avoid shearing in the perturbation, we should impose Re(w,) 
to have the same constant value in all layers and just look at the 
variation of the growth rate, Re(w,) is related to the pattern speed 
of the perturbation. As a result, we would expect that the growth 
rate depends on radius and is maximal where the shear is largest. 
The column will then split into weakly (almost no shear in the 
flow) and strongly (large velocity gradients) unstable regions. 



3.5. Algorithm 

The eigenvalue problem, Eq. ( |52] ). is solved by standard numer- 
ical techniques. We have implemented a shooting method as de- 
scribed in Petri (2007a) by replacing the eigenvalue Eq. (54) 
there, by the new ordinary differential equation, Eq. ( |52l ) and us- 
ing the jump conditions Eq. ( |56] | and ( |57] |. when inner and outer 
walls are present. 

For the pulsar electrosphere, the situation is very similar, ex- 
cept that no calculation is performed in region III. The bound- 
ary condition for outgoing waves is applied at R2, see Eq. 
Actually for pulsars, we compare both boundary conditions. 



4. Algorithm check 

In order to check our algorithm in different configurations, 
we compute the eigenvalues for both, a low-density relativistic 
plasma column and in the limiting case of a relativistic planar 
diode geometry. For some special density profiles, the exact ana- 
lytical dispersion relations are known and used as a starting point 
or for comparison with the numerical results obtained by our al- 
gorithm. 

4.1. Low-density plasma column 

In cylindrical geometry, an exact analytical solution of the 
dispersion relation can be found in the l ow-density and non- 
relativistic limit, the diocotron instability (IDavidsonlll990l) . We 
use these results to check our algorithm in cylindrical coordi- 
nates, as was already done in Petri (2007a). 

For completeness, we briefly recall the main characteristics 
of the configuration. The magnetic field is constant and uniform 
outside the plasma column, namely - Bq, for r > R2. We 
solve self-consistently the Maxwell equations in the space be- 
tween the inner and the outer wall, Ri < r < R2. The parti- 
cle number density and charge density are constant in the whole 
plasma column such that 



f , Wi<r<Ri 

Pt(r) = I po = no ? = const , Ri <r <R2 
, R2<r<W2 



(63) 



We first check our fully relativistic and electromagnetic code, 
which includes particle inertia, in the non-relativistic and low- 
density limit. By non-relativistic, we mean a maximum speed at 
the outer edge of the plasma column of a: 10"^ and a density such 
that the self-field parameter is weak, Se - 10"^. When numerical 
values of the self-field parameter are given, we mean the value it 
takes at the outer edge of the plasma column. Thus ig should be 
understood as s^^ir - R2). In this limit, the electric drift approxi- 
mation is excellent and the non-relativistic diocotron regime ap- 
plies. The eigenvalues are therefore given by Eq. (68) of lP6tril 
(2007a). 

We computed the growth rates for several geometrical con- 
figurations of the plasma column, varying R12 and Wi,2. Some 
typical examples are presented in Table [T] To summarise, the 
relative error in the real and imaginary part of the eigenvalues 
compared to the exact analytical solution is 



Sim - 



Re (OJ) - Re (Wexact) 



Re (Wexact) 
Im (CJ) - Im (Wexact) 



Im(Wexact) 



< 10" 



< 10" 



(64) 
(65) 
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iVHJU-C L 


rl, 




^num 




^Im 


2 


0.4 


0.5 


1.8865e-03 + 3.5878e-04 i 


6.5774e-08 


4.1221e-08 


3 


0.4 


0.5 


2.7284e-03 + 1.1336e-03 i 


7.8227e-08 


1.3019e-08 


4 


0.4 


0.5 


3.608 le-03 + 1.4944e-03 i 


8.6291e-08 


1.2482e-08 


5 


0.45 


0.5 


2.3766e-03 + 1.3515e-03 i 


1.9263e-08 


3.1651e-09 


7 


0.45 


0.5 


3.325 le-03 + 1.7069e-03 i 


2.2522e-08 


4.2726e-09 



Table 1. Numerical eigenvalues Wnum and relative errors for the low-density and non-relativistic plasma column, for different modes I 
and different aspect ratios, d\ - R1IW2, and di - RilWi. 



The precision is excellent, reaching 7 digits at least. Our algo- 
rithm computes quickly and accurately the eigenvalues in cylin- 
drical geometry with vacuum gaps between the plasma column 
and the walls. The eigenvalues obtained in this example are good 
initial guesses to study the relativistic problem in the low speed 
limit. 



As in iPetril (I2007al) . the influence of the relativistic and 
electromagnetic effects are investigated by slowly increasing 
the maximal speed at the outer edge of the plasma column, 
jS2 = Ri^iRiilc. We also allow for outgoing electromagnetic 
wave radiation at the outer edge of the column. 

Two cases are presented in Fig.[T] The first one, Fig.[Tt, has 
w = W1/W2 = 0.1, di = R1IW2 = 0.4 and d2 = R2IW2 = 0.5 
whereas the second one. Fig. [TJi, has w - 0.1, d\ - 0.45 and 
d2 - 0.5. For non-relativistic speeds, ySa 1, the eigenvalues of 
Sect. 14. II are recovered. The thinner the plasma layer, the larger 
the number of unstable modes, respectively 5 and 12 unstable 
modes. In both cases, the growth rate starts to be altered when- 
ever J02 > 0.1. In any case, for very high speeds, ySa ~ 1, all the 
modes stabilise because the growth rate vanishes. 

Removing the outer wall, and replacing it by outgoing wave 
boundary conditions does not significantly affect the instability. 
Except for the mode I - 2, the changes become perceptible only 
for relativistic speeds, yS2 > 0.7. 

Note that th ese results slightly differ from those presented 
in iPetril (12007 a') because in the former case the density is kept 
constant while in the latter case the diocotron frequency, (Wd, is 
kept constant, both being related by the Lorentz factor of the 
flow. 

A simplistic approach to understand the stabilisation process 
is as follows. In the low-density limit, the magnetron instability 
reduces to the diocotron regime. Assuming a constant charge 
density within the plasma, unstable oscillations are generated 
when the surface waves at both edges of the plasma column 
interact in a constructive way, i.e. that oscillation have to syn- 
chronise. For high velocities within the fluid, the relative phase 
speed between the two waves increase, mutual synchronisation 
beco mes difficult to r each, reducing the growth rate of the insta- 
bility (|^M^[i966|). 

The qualitative physical nature of the diocotron/magnetron 
instability can be expressed via t he m ore famil iar Kelvin- 
Helmholtz or two-stream instability ( MacF arlane & Havi Il950t 
iTrivelpiece & Gouidlll959l : lBuneman et aUil96 6). For instance, 
the tw o-stream instability is investigated in 'Delcroix & Beri 
(II994I) . It is found (in Cartesian geometry) that while in the long 
wavelength limit, the growth rate is proportional to the velocity 
difference between the two beams, they are reduced for suffi- 
ciently large velocity shear until the instability disappears. This 
stabilisation mechanism applies also to the diocotron/magnetron 
instability and can already occur at modest speeds, much less 
than the speed of light c. 



4.2. Influence of curvature 

The influence of the curvature is also studied by taking the limit 
of the planar diode geometry. The curvature of the plasma col- 
umn is then increased to investigate the evolution of the growth 
rates. 

The effective aspect ratio of the plasma layer is conveniently 
described by the parameter 



R\ 



R2 - Rl 



(66) 



Using the equilibrium electric field profile indicated in Sect. 12.51 
in the limit of small curvature corresponding to large aspect ra- 
tio, A -Hoo, the eigenvalue problem and equilibrium configu- 
ration is described by the relativistic planar diode. 

We show the evolution of the growth rate in the non- 
relativistic limit - lO""*, Fig. |2^, and in the relativistic case 
/32 ~ 0.2, Fig.|2]3. The aspect ratio has a drastic influence on the 
growth rate. For large values of A » 1, all unstable modes are 
stabilised, in both non-relativistic and r elativistic flows. These 
results agree with those found in IPetril ( 12007 a ). Note that the 
drift speed has only a negligible effect on the growth rate com- 
pared to the aspect ratio. 

4.3. Electrosphere 

Finally, we checked the results obtained for the electrosphere, 
with outer wall or outgoing waves. The rotation profile is chosen 
to mimic the rotation curve obtained in the 3D electrosphere. To 
study the influence of the relativistic ef fects, we take the same 
profiles as those given in 'Petri' ('2007a"F). We remind that dif- 
ferent analytical expressions for the radial dependence of Q are 
chosen by mainly varying the gradient in differential shear as 
follows 



Q.(r) = (2 + tanh[a' (r - r,,)] e"^'*) 



(67) 



Q., is the neutron star spin and r is normalised to the neutron star 
radius, R,. The values used are listed in Table|2] The angular ve- 



Q. 


a 




'0 




1.0 


5 X 10-' 


6.0 




0.3 


5 X 10-5 


10.0 



Table 2. Parameters for the rotation profiles used to mimics the 
azimuthal velocity of the plasma in the electrospheric disk. 



locity, and the corresponding particle density number, are shown 
respectively in Fig.[3t and[3]3. In both cases, Q starts from coro- 
tation with the star Q = 0» followed by a sharp increase around 
r - 6 for D.2 and a less pronounced gradient around r = 10 for 
Q3. Finally the rotation rate asymptotes twice the neutron star 
rotation speed for large radii. 
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Fig. 1. Stabilisation of the magnetron instability in the low-density limit. We retrieve the results of the diocotron instability. The 
plasma density is constant in the column, po - est. The geometric aspect ratios are, for Fig. (a), w - 0.1, d\ - 0.4 and d2 - 0.5 
and, for Fig.(b), w - 0.1, d\ - 0.45 and d2 - 0.5. The growth rates are normalised to the diocotron frequency, wd- The curves 
corresponding to the outgoing wave boundary conditions contain twice as much symbols as the conducting wall boundary conditions 
in order to distinguish them. 
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Fig. 2. Effect of the cylindrical geometry on the growth rate of the instability. The growth rate are normalised to the rotation at the 
outer edge of the plasma column Q(/?2) and plotted versus the aspect ratio A, in logarithmic scale. In Fig. (a), the speed at the outer 
edge of the column is j32 - lO^^* whereas in Fig.(b), it is P2 ~ 0.2. 
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Fig. 3. Two choices of differential rotation curves in the plasma column for the cylindrical pulsar electrosphere, Q.2 with strong shear 
(in green) and Q3 with weak shear (in blue). Fig. (a) and the corresponding charge density, Fig.(b). 



The computed growth rates, normalised to the speed of the cylindrical coordinates, for various equilibrium density profiles, 
neutron star, are shown in Fig. HI for the profile D.2, Fig.H^, and electric field, and velocity profiles. AppUcation to pulsar's elec- 
for fij, Fig.|4j3. trosphere is also discussed. 



5. RESULTS 

We demonstrated that our numerical algorithm gives accurate 
results in the non-relativistic cylindrical geometry for which 
we know an analytical expression of the eigenfrequencies. 
Moreover, the transition to the relativis tic regime gives the same 
results as those shown in lPetril(l2007al) . In this section, we com- 
pute the eigenspectra of the relativistic magnetron instability in 



5.1. Influence of the self -field 

We study the influence of the self-electromagnetic field on the 
instability. When the self-field is negligible, i.e. in the low- 
density limit, ie « 1 , the instability is well approximated by the 
diocotron regime, which means in the electric drift approxima- 
tion. However, when the density of the plasma increases, ie ~ 1, 
it induces a strong electric field which would lead to superlumi- 
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Fig. 4. Evolution of the growth rates, lm(a>), normalised to the spin of the neutron star Q,, for increasing maximal speed of the 
column, j62. For the profile Q.2, in Fig. (a), outer conducting wall are compared with outgoing wave boundary conditions. For the 
latter, the curves contain twice as much symbols as for the former boundary conditions. Note that for the mode I - I, the curves 
overlap and cannot be distinguished. For the profile Q3, in Fig.(b), only the wall boundary conditions are plotted. 



nal motion in the drift approximation. Therefore, particle inertia 
comes into play and imposes a motion which departs signifi- 
cantly from the electric drift and modifies the behavior of the 
instability. 

We considered the plasma in cylindrical geometry, confined 
by some external experimental electromagnetic device between 
the inner and outer wall. The external applied magnetic field and 
the density profile are specified as initial data. Thus, we deduce 
the velocity profile by solving the non-linear Volterra integral 
equation, Eq. ( |26] |, as explained in Sect. 12.41 

We study the influence of particle inertia in both, the non- 
relativistic and relativistic regime. As a starting point, we take 
the diocotron regime, «: 1, and slowly increase the influence 
of the self-field parameter, Sg. 

Examples of growth rates are shown in Fig. |5] A non- 
relativistic flow (J32 ~ 10"^) is shown in Fig.|5^ whereas a mildly 
relativistic flow (J32 ~ 0.19) is shown in Fig.jSj^- In both cases, 
increasing the self-field stabilises the magnetron regime. Close 
to the maximum value of the self-field parameter, the Brillouin 
zone above which no equilibrium configuration exists because 
the defocusing electric field is to strong compared to the mag- 
netic focusing field, the instability has almost completely disap- 
peared. 

For completeness, the evolution of the shape of the eigen- 
functions of the mode / = 5 corresponding to the eigenvalues 
in Fig. |5h is shown in Fig. |6] When the stabilisation process 
becomes important, the real and imaginary parts of the eigen- 
function possess significant negative values compared to the low 
self-field parameter case. 

5.2. Electrosphere 

The electrospheric non-neutral plasma is confined by the rotat- 
ing magnetised neutron star. The most important feature is the 
velocity profile in the plasma column. For simplicity, here, we 
assume that no vacuum gaps exist between the plasma and the 
inner wall, Wi - Ri. The inner wall depicts the perfectly con- 
ducting neutron star interior. The plasma is in contact with the 
neutron star surface. The outer wall can be suppressed depending 
on the outer boundary conditions. For instance, we expect elec- 
tromagnetic wave generation, leading to a net outgoing Poynting 
flux carrying energy to infinity. Both situa tions will be consid- 
ered. We generalise the study presented in iPetril (l2007al) by in- 
cluding particle inertia. 

To remain fully self-consistent, we only consider an uniform 
applied external magnetic field. The growth rates for the rotation 



curve Q2 for each mode I are shown in Fig. |7h and those for Q3 
in Fig.lTj). 

For Q.2, starting from the diocotron regime, s^. I, the 
instability is first enhanced by the particle inertia effect, un- 
til ie ~ 0.25. Then, the growth rates begin to decrease until 
they vanish for the maximum allowed self-field parameter, cor- 
responding to the last existing equilibrium configuration. 

For Q3, no increase is observed, the stabilisation eff'ect pro- 
ceeds as soon as increases. 

Finally, we compute the evolution of the spectrum of the 
magnetron instability when going to the relativistic regime for 
a significant initial self-field parameter, x 0.16. We start with 
a non-relativistic rotation profile such that /32 I and slowly in- 
crease R2 (as well as Wi,W2,Ri to maintain their ratio constant) 
in order to approach the speed of light for the maximal rotation 
rate of the plasma column. 

Results for the growth rates of the velocity profiles Q.2 and 
Q3 are shown in Fig. [8^ and b respectively. 

For completeness, the evolution of the shape of the eigen- 
functions of the mode / = 5 corresponding to the eigenvalues 
in Fig. [8^ is shown in Fig. |9] When the stabilisation process 
starts, the real and imaginary parts of the eigenfunction oscil- 
late. Moreover, while these functions remain positive for small 
speeds, they change sign whenever the diminishing of the growth 
rates become significant. This behavior was already observed in 
the study of the influence of the self-field parameter. Sect. 15.11 

6. CONCLUSION 

We developed a numerical code to compute the eigenspectra and 
eigenfunctions of the magnetron instability in cylindrical geom- 
etry, including particle inertia, as well as electromagnetic and 
relativistic effects, in a fully self-consistent way. It is thus possi- 
ble to study the behavior of the plasma in the vicinity of the light 
cylinder, i.e. where the particle kinetic energy becomes compa- 
rable to the magnetic field energy density. 

Unstable modes are computed for a uniform external ap- 
plied magnetic field and arbitrary velocity, density and electric 
field profiles. The resulting equilibrium configuration is com- 
puted self-consistently, according to the cold-fluid and Maxwell 
equations. Application of the code to a plasma column as well 
as to the pulsar electrosphere have been shown. In both cases, 
the magnetron regime gives rise to instabilities that become less 
and less unstable when the self-field parameter of the flow, ig, 
approaches its maximum allowed value, < ie.max- Whereas 
the growth rates can be comparable to the rotation period of 
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Fig. 5. Influence of the self-field parameter, s^^, on the magnetron instability, for a non-relativistic flow (J52 ~ 10 ^) in Fig. (a) and a 
relativistic flow (J32 ~ 0.19) in Fig.(b). The geometric aspect ratios are, w - Q.l,di = 0.4, and d2 = 0.5, the same as those in Fig.[T^. 
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Fig. 6. Real and imaginary part of the eigenfunction, respectively Fig. a and b, for the mode / = 5 in Fig.|5^. The real part of the 
eigenfunction is normalised such that the maximum in absolute value is unity. The value of the self-field parameter for each 
eigenfunction is shown in legend. 
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Fig. 7. Evolution of the growth rates, Im(w), of the profiles Q2, Fig. (a), and Q3, Fig.(b), for increasing self-field parameter, i.e. 
transition from the diocotron to the magnetron regime. 
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Fig. 8. Evolution of the growth rates Im(w) of the profiles Q2, Fig. (a), and Q3, Fig.(b), for a significant initial self-field parameter, 
Se ~ 0.16. 
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Fig. 9. Real and imaginary part of the eigenfunction, respectively Fig. a and b, for the mode / = 5 in Fig. [8^. The real part of 
the eigenfunction is normalised such that the maximum in absolute value is unity. The value of the maximal speed [32 for each 
eigenfunction is shown in legend. 



the neutron star in the non-relativistic limit, it is found that for 
special rotation profiles, the magnetron instability is completely 
suppressed in the relativistic regime, as was already the case in 
the diocotron limit. 

However, the magnetron instability is not restricted to the 
electrosphere. For instance, it is believed that the pulsed radio 
emission emanates from the magnetic poles close to the neutron 
star surface. Highly relativistic and dense leptonic flows are ac- 
celerated at these polar caps, forming non-neutral plasma beams 
able to radiate coherent electromagnetic waves by several pro- 
cesses like the cyclotron maser or the free electron laser insta- 
biUty. We will address this problem as well as the eflfect of finite 
temperature in a forthcoming paper. 
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